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ABSTRACT 

\ We study the distribution of dark matter in the Draco dwarf spheroidal galaxy by 

modelling the moments of the line-of-sight velocity distribution of stars obtained from 
new velocity data of Wilkinson et al. The luminosity distribution is approximated by 
a Sersic profile fitted to the data by Odenkirchen et al. We assume that the dark mat- 
ter density profile is given by a formula with an inner cusp and an outer exponential 
cut-off, as recently proposed by Kazantzidis et al. as a result of simulations of tidal 
stripping of dwarfs by the potential of the Milky Way. The dark matter distribution 
is characterized by the total dark mass and the cut-off radius. The models have arbi- 
trary velocity anisotropy parameter assumed to be constant with radius. We estimate 
the three parameters by fitting both the line-of-sight velocity dispersion and kurtosis 
profiles, which allows us to break the degeneracy between the mass distribution and 
velocity anisotropy. The results of the fitting procedure turn out to be very different 
depending on the stellar sample considered, that is on our choice of stars with dis- 
crepant velocities to be discarded as interlopers. For our most reliable sample, the 
model parameters remain weakly constrained, but the robust result is the preference 
for weakly tangential stellar orbits and high mass-to-light ratios. The best-fitting total 
mass is then 7 x 10 7 Mq, much lower than recent estimates, while the mass-to-light 
ratio is M/Ly = 300 and almost constant with radius. If the binary fraction in the 
stellar population of Draco turns out to be significant, the kurtosis of the global veloc- 
ity distribution will be smaller and the orbits inferred will be more tangential, while 

the resulting mass estimate lower. 
• i— i . 
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1 INTRODUCTION 

Dwarf spheroidal galaxies are currently believed to be among 
the most dark matter dominated objects in the Universe 
(with mass-to-light ratios up to few hundred solar units) 
and therefore are of immediate interest when trying to test 
present theoretical predictions concerning the dark matter 
profiles. As dwarfs, they are also important for theories of 
structure formation, which suffer from the overabundance of 
subhaloes (Klypin et al. 1999). 

The Draco dwarf spheroidal galaxy (hereafter, Draco) 
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is a generic example of this class of objects in the neigh- 
bourhood of the Galaxy and has been already the subject of 
extensive study. Recently, new radial-velocity observations 
of Draco have yielded a larger kinematic sample (Kleyna et 
al. 2002; Wilkinson et al. 2004) and the Sloan Digital Sky 
Survey has yielded a better determination of the luminosity 
profile and shape (Odenkirchen et al. 2001). Draco has also 
been examined for the presence of tidal tails which could 
indicate a kinematically perturbed state. For a long time it 
has been debated (e.g. Klessen & Kroupa 1998; Klessen & 
Zhao 2002) whether the large velocity dispersion measured 
in Draco is a result of significant dark matter content or a 
signature of tidal disruption by the Milky Way. New analy- 
ses by Piatek et al. (2002) and Klessen, Grebel & Harbeck 
(2003) show that there is no evidence that tidal forces from 
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the Milky Way have disturbed the inner structure of Draco. 
In particular, Klessen et al. (2003) argue that models with- 
out any dark matter are unable to reproduce the narrow 
observed horizontal branch width of Draco. As we will show 
below, the truth may lie between these two pictures: Draco 
most probably is strongly dark matter dominated but still 
its outer edges could be affected by tidal interaction with 
the Milky Way producing a population of stars unbound to 
the dwarf. 

One of the major methods to determine the dark matter 
content of galaxies and clusters is to study the kinematics 
of their observable discrete component: stars or galaxies, 
respectively. The method is based on the assumption that 
the positions and velocities of the component members trace 
the Newtonian gravitational potential of the object. If the 
object is in virial equilibrium the relation between those is 
described by the Jeans formalism. In the classical approach, 
such an analysis is restricted to solving the lowest order 
Jeans equation and modelling only the velocity dispersion 
profile. In previous papers (Lokas 2002; Lokas & Mamon 
2003), building on the earlier work by Merrifield & Kent 
(1990) and van der Marel et al. (2000), we have extended 
the formalism to include the solutions of the higher-order 
Jeans equations describing the fourth velocity moment, the 
kurtosis. 

The formalism has been successfully applied to study 
the dark matter distribution in the Coma cluster of galax- 
ies (Lokas & Mamon 2003). We have shown that, for a re- 
stricted class of dark matter distributions motivated by the 
results of cosmological iV-body simulations, the joint analy- 
sis of velocity dispersion and kurtosis allows us to break the 
usual degeneracy between the mass distribution and velocity 
anisotropy and constrain the parameters of the dark mat- 
ter profile. Recently we have tested the reliability of this 
approach against a series of iV-body simulations (Sanchis, 
Lokas & Mamon 2004) and verified that the method allows 
for typically accurate determinations of the mass, concen- 
tration and velocity anisotropy (albeit with a large scatter 
in concentration). In the present paper, encouraged by these 
successes, we apply the method to constrain the dark matter 
distribution in the Draco dwarf spheroidal galaxy. 

The paper is organized as follows. In Section 2 we sum- 
marize the Jeans formalism for the calculation of the mo- 
ments of the line-of-sight velocity distribution. In Section 3 
we estimate the moments from the velocity data and discuss 
the associated difficulties caused by the possible presence of 
interlopers and binary stars. Section 4 describes our models 
for the luminosity and dark matter distributions. Results of 
fitting the models to the data are presented and discussed in 
Section 5. In Section 6 we comment on the possible origin of 
unbound stars in Draco and the concluding remarks follow 
in Section 7. 



2 VELOCITY MOMENTS FROM THE JEANS 
FORMALISM 

The Jeans formalism (e.g. Binney & Tremaine 1987) relates 
the velocity moments of a gravitationally bound object to 
the underlying mass distribution. We summarize here the 
formalism, as developed in Lokas (2002) and Lokas & Ma- 
mon (2003). The second of and fourth-order vf radial ve- 



locity moments obey the Jeans equations 
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where v is the 3D density distribution of the tracer popula- 
tion and $ is the gravitational potential. The second equa- 
tion was derived assuming the distribution function of the 
form f(E,L) — fo(E)L~ 2/3 and the anisotropy parameter 
relating the angular a$ and radial oy velocity dispersions 
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to be constant with radius. We will consider here — oo < /3 < 
1 which covers all interesting possibilities from radial orbits 
(/3 = 1) to isotropy (/3 = 0) and circular orbits (/3 — > — oo). 
The solutions to equations are 
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Projecting along the line of sight we obtain the observable 
quantities 
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where 

g(r, R,0) = l- 2/3— + — 

I{R) is the surface distribution of the tracer and 7? is the 
projected radius. 

Introducing equations 0-10 into ©-0 an d inverting 
the order of integration, the calculations of <ti os and vf os can 
be reduced to one- and two-dimensional numerical integrals 
respectively (see the appendix of Mamon & Lokas 2005 for 
a simpler expression for <ti os ). In the following we will refer 
to the fourth moment in the form of kurtosis 
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whose value is 3 for a Gaussian distribution. 



3 VELOCITY MOMENTS FROM 
OBSERVATIONS 

Fig. shows the heliocentric velocities and projected dis- 
tances from the centre of 207 stars with good velocity 
measurements from Wilkinson et al. (2004). The distances 
are calculated assuming that the centre of Draco is at 
RA=17 h 20 m 13.2 s , Dec=57°54'54" (J2000) (Odenkirchen et 
al. 2001). The sample was obtained from the original sample 
of 416 velocities stars observed in the direction of Draco by 
cutting out all objects with velocities differing by more than 
39 km s _1 from Draco's mean velocity of -291 km s" 1 . This 
choice was motivated by the assumption that the velocity 
distribution is Gaussian with a maximum dispersion of 13 
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Figure 1. Heliocentric velocities versus distances from the centre 
of Draco for 207 stars. The solid (dashed) curves separate 203 
(189) stars counted as members of the galaxy from the 4 (18) 
supposed interlopers. 



km s _1 . We should however still take into account the possi- 
bility that some of the 207 stars are interlopers unbound to 
the gravitational potential of the galaxy. Prada et al. (2003) 
discussed the issue thoroughly in the case of isolated galax- 
ies in the Sloan Digital Sky Survey (SDSS) and showed that 
the presence of interlopers can affect significantly the veloc- 
ity dispersion and therefore the inferred mass distribution 
in the galaxy. This effect is especially important in the case 
of small samples like the present one. 

The problem has been also discussed by Mayer et al. 
(2001) who showed (see Fig. 26 of the published version of 
their paper) how the velocity dispersion profiles of the dwarf 
spheroidal galaxies can be affected if the line of sight is along 
one of the tidal tails: while the intrinsic dispersion profile is 
declining, the observed one, due to the presence of unbound 
stars in the tail, is much higher and increasing. They state 
explicitly that a velocity dispersion that increases outside 
the core of the dwarf should be considered as a result of 
possible contamination by tidal tails. We address the issue 
in more detail in Section 6. 

Judging by the shape of the diagram shown in Fig. Q 
and comparing it with the expected appearance of gravita- 
tionally bound object in velocity space we immediately see 
that there are at least four stars with discrepant velocities: 
they have been separated from the main body of the galaxy 
by the solid lines. An even more stringent, but still arbitrary 
choice, with 14 more stars excluded, is shown with dashed 
lines. The lines were drawn symmetrically with respect to 
the mean systemic velocity of —291 km s , as estimated 
from the total sample with 207 stars. Since we have no possi- 
bility of determining a priori which stars are actually bound 
to the galaxy we calculate the velocity moments using the 
three samples with 207, 203 and 189 velocities obtained in 
this way. 

In the estimation of the moments, we follow the Ap- 
pendix in Lokas & Mamon (2003) neglecting the small errors 



in the velocity measurements. We divide the sample into 5 
radial bins each containing 37-43 stars (depending on the 
sample: 207 stars are divided into bins with 4 x 41 + 43 ve- 
locities, 203 stars into 4 x 41 + 39 velocities and 189 stars 
into 4 x 38 + 37 velocities). The most natural estimators of 
the variance and kurtosis from a sample of n line-of-sight 
velocity measurements Vi are 
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(12) 



is the mean of stellar velocities in the sample. The distribu- 
tions of these estimators for our binning of stars, i.e. when 
n ~ 40, can be investigated by running Monte Carlo sim- 
ulations selecting AT = 10 4 times n = 40 numbers from 
a Gaussian distribution with zero mean and dispersion of 
unity (see Lokas & Mamon 2003). One can then construct 
unbiased and nearly Gaussian-distributed estimators of line- 
of-sight velocity dispersion s and kurtosis-like variable k 



k = 



(— 

\n - 1 



K 



1/10 



(13) 
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The factor n — 1 in equation il l Ml is the well known cor- 
rection for bias when estimating the sample variance, valid 
independently of the underlying distribution. In 1141 the fac- 
tor 3/2.75 corrects for the bias in the kurtosis estimate, i.e. 
unbiased estimate of kurtosis is K' — 3if/2.75, while the 
rather complicated function of K' assures that the sampling 
distribution of k is approximately Gaussian. We find that 
the standard errors in the case of s are of the order of 11 
percent while in the case of k are approximately 2 percent. 
In the following we assign these errors to our data points. 

We checked that even for weakly non-Gaussian veloc- 
ity distributions, such as following from the present data, 
estimators s and k are Gaussian-distributed to a very good 
approximation and very weakly correlated with the corre- 
lation coefficient \g\ < 0.07. We can therefore assume that, 
to a good approximation, all of our data points measuring 
velocity dispersion and kurtosis will be independent, which 
justifies the use of standard \ 2 minimization to fit the mod- 
els to the data. 

Apart from the sampling distributions discussed above, 
our estimates of the velocity moments can be affected by the 
presence of binary stars. Contrary to the studies of galaxy 
clusters, where galaxy pairs are rare and can be rather eas- 
ily identified and eliminated from the sample, in the case 
of stellar systems, such as Draco, binaries can be identified 
only through long-term monitoring of the whole stellar sam- 
ple. In the sample of Wilkinson et al. (2004) such repeated 
measurements exist only for 66 stars and among those only 
a few are probable binaries, i.e. the differences between their 
velocity measurements are too big to be accounted for by the 
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Figure 2. Velocity moments calculated from the three samples 
of 207 (solid lines), 203 (dashed lines) and 189 stars (dotted lines) 
and different binary fraction 7 = 0, 0.25 and 0.5. In each triplet 
the highest line corresponds to 7 = (no correction for binaries). 



measurement errors. Given the uncertainties in our knowl- 
edge of the binary population in Draco, we assume that the 
binary fraction is negligible (7 = 0), but we discuss briefly 
below how a significant amount of binaries could affect the 
results. 

De Rijcke & Dejonghe (2002) studied the influence of 
the binary population on the observed line-of-sight velocity 
moments and found (see their Fig. 15) that for stellar sys- 
tems with <Ti os ~ 9 km s -1 , such as Draco, the observed 
(Gaussian) dispersion is mildly affected by the presence of 
binaries in comparison with the intrinsic stellar velocity dis- 
persion. However, the kurtosis value turns out to be more 
affected. Inverting the formulae (55) in De Rijcke & De- 
jonghe (2002) we find that given the line-of-sight velocity 
moments of the binary population, <7b and Kb we can calcu- 
late the intrinsic moments of the stars cri and Ki from the 
observed ones o a and k 



<7j 



2 

7°b 



7^ 



67 



(15) 
(16) 



We adopt the values of <Tb = 2.87 km s" 1 and k\, = 86.89, as 
obtained by De Rijcke & Dejonghe (2002) for their standard 
binary population model. Having calculated the observed 
moments a = s and k q = K' — 3K/2.75 from equations 
l|10|l - l|13| l we obtain the intrinsic values from equations l|15|l - 

The uncertainties in the determination of the velocity 
moments are illustrated in Fig. [2] The solid lines join the 
moments calculated in five bins with all 207 velocities, the 
dashed and dotted ones with 203 and 189 velocities respec- 
tively. The three curves of each kind correspond to the bi- 
nary fraction 7 = 0, 0.25 and 0.5, from the highest to the 
lowest, i.e. correcting for binaries decreases the values of 
moments. The estimates overlap for the first three bins of 
samples with 207 and 203 stars because none of the supposed 



interlopers fall there, the small differences in the case of the 
sample with 189 stars are due to different number of stars 
per bin. For the two outer bins the situation is however very 
different, with all 207 stars the velocity dispersion is much 
larger than for 203 and 189 stars. We note that although 
the profile for 207 stars decreases in the last bin, there is no 
sudden drop in velocity dispersion as found by Wilkinson et 
al. (2004). This is mainly due to our larger number of stars 
per bin. 

As demonstrated by Fig. [5] except for the most numer- 
ous sample with 207 stars, the kurtosis values typically fall 
below the Gaussian value of 3. Such kurtosis signifies a cen- 
trally flattened velocity distribution in comparison with the 
Gaussian. It is worth noting that the behaviour of the kur- 
tosis profile for the sample with 189 stars (close to Gaussian 
value in the centre and lower values at larger distances) is 
similar to the one found in the case of distribution of galaxy 
velocities in the Coma cluster (Lokas & Mamon 2003). De- 
creasing kurtosis values seem also to be typical properties of 
the projected velocity distributions of particles in simulated 
dark matter haloes as demonstrated by the studies of San- 
chis et al. (2004), Kazantzidis, Magorrian & Moore (2004) 
and Diemand, Moore & Stadel (2004). 

Our purpose now is to reproduce the observed profiles 
of velocity moments using models described in the next Sec- 
tion. 



4 DISTRIBUTIONS OF STARS AND DARK 
MATTER 

4.1 Stars 

The distribution of stars is modelled in the same way as in 
Lokas (2001, 2002) i.e. by the Sersic profile (Sersic 1968) 



I(R) = 7 exp[-( J R/ J R s ) 



I/mi 



(17) 



where Iq is the central surface brightness and Rs is the char- 
acteristic projected radius of the Sersic profile. The 3D lu- 
minosity density v(r) is obtained from I(R) by deprojection 



dl dR 



dR^/R2~T 



(18) 



For a wide range of values 1/2 < m < 10 found to fit 
the light distribution in elliptical and spheroidal galaxies an 
excellent approximation for v(r) is provided by (Lima Neto, 
Gerbal & Marquez 1999) 



/ r(2m) 



exp 



/ r y/,-, 



(19) 



2R s r[{3-p)m] 
p = 1.0 - 0.6097/m + 0.05463/m 2 . 

The mass distribution of stars following from I19|l is 

where Ms = TL to t is the total mass of stars, with T=const 
the mass-to-light ratio for stars, Ltot the total luminosity 
of the galaxy, and j(a,x) = J Q e~ t t a ~ 1 dt is the incomplete 
gamma function. 
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Table 1. Adopted parameters of the Draco dwarf. 



parameter 


value 


luminosity Ltot.v 


2.2 X 1O 5 L 


distance d 


80 kpc 


distance modulus m — M 


19.5 


stellar mass-to-light ratio 


3M Q /L Q 


stellar mass M$ 


6.6 x 10 s M Q 


Sersic radius Rg 


7.3 arcmin 


Sersic parameter m 


0.83 



Odenkirchen et al. (2001) find that their data for the 
surface density of stars in Draco are best fitted (for sam- 
ple SI) with 1/m = 1.2 (m = 0.83) and the Sersic radius 
Rs — 7.3 arcmin. We will adopt these values here assum- 
ing that the population of stars is uniform and therefore the 
distribution of stars can be translated to the surface lumi- 
nosity in any band just by adjusting the normalization 7o. 
We note that the surface density distribution estimated by 
Odenkirchen et al. (2001) agrees well with the independent 
estimates by Piatek et al. (2002) from a different data set. In 
particular, both groups find the distribution to fall sharply 
at large distances thus contradicting the earlier estimates of 
Irwin & Hatzidimitriou (1995) who found a flattening profile 
which could point towards the presence of tidal tails. 

We adopt the apparent magnitude of Draco in the i- 
band m% = 10.49 following Odenkirchen et al. (2001). Us- 
ing their formula (8) relating i and / magnitudes, and as- 
suming V — I — 1 from the colour-magnitude diagram in 
Fig. 5 of Bonanos et al. (2004) we find i — I = 0.52 in 
agreement with Table 3 of Fukugita, Shimasaku & Ichikawa 
(1995) and V — i — 0.48. Taking a distance estimate for 
Draco of d = 80 kpc (Aparicio, Carrera & Martinez-Delgado 
2001), corresponding to the distance modulus m—M — 19.5, 
we obtain the total luminosity of Draco in the U-band of 
itot.v = 2.2x 10 5 Lq, which is somewhat larger than the ear- 
lier estimate of Ltot,V = 1.8 x 10 5 Lq by Irwin & Hatzidim- 
itriou (1995). We note that the adopted distance agrees well 
with the most recent estimate by Bonanos et al. (2004) who 
found m — M = 19.4 once an appropriate (for their method 
of distance scale calibration) correction of 0.1 mag is applied. 

In order to estimate the mass in stars, we need to adopt 
the stellar mass-to-light ratio in l/-band, which for dwarf 
spheroidals is believed to be similar to that of globular clus- 
ters Tv = (1 — 3)Mq/Lq. Although the age and metallic- 
ity of Draco stellar population are reasonably well known, 
the predicted values of Tv are still uncertain and different 
stellar population synthesis codes and different initial mass 
functions tend to produce results discrepant at least by a fac- 
tor of 2 (Schulz et al. 2002). Since we do not want to overesti- 
mate the amount of dark matter in our modelling, and keep- 
ing in mind that there is probably a significant contribution 
from ionized gas in Draco (up to 2 x 10 5 Mq, Gallagher et 
al. 2003) we take a higher value Tv = 3Mq/Lq so that our 
total stellar mass will be Ms = TyLtot.v = 6.6 x 10 Mq. 
All the parameters of Draco discussed in this Section are 
summarized in Tabled 

We also note that the ellipticity of the surface density 
distribution of stars in Draco is low, e = 0.3 (Odenkirchen 
et al. 2001; Piatek et al. 2002), so the galaxy is likely to 
be close to spherical, which makes the use of our spherical 



models justified and reduces the chance that non-sphericity 
will affect our results (see Sanchis et al. 2004). 

4.2 Dark matter 

Recently, Kazantzidis et al. (2004) performed a series of N- 
body simulations of the evolution of a subhalo in the field 
of a larger halo of a galaxy. Their Draco-like object ini- 
tially had an NFW density distribution (Navarro, Frenk & 
White 1997) but due to tidal interactions with a galaxy sim- 
ilar to the Milky Way lost some of its mass. The resulting, 
tidally stripped, density profile retained the initial r~ cusp 
at the centre but developed an exponential cut-off at large 
distances and could be well approximated by a formula 



p d (r)=Cr l exp^-^ 



(21) 



where r b is the break radius at which the cut-off occurs. The 
dark mass distribution associated with the density profile 
ED is 



M d (r) = M D 



(22) 



1 _ exp 1 + - 

L V r h J \ r b/ 

which we have normalized with Mb, the total dark mass of 
the halo. (Note that contrary to the standard NFW dis- 
tribution, for this profile the mass converges.) With this 
normalization, the constant in equation 12 U becomes C = 
M D /(47rrg). 

It is still debated whether tidal interaction of dwarfs 
with gravitational potential of their host galaxies is likely to 
flatten the inner profile of their halo. Although Kazantzidis 
et al. (2004) find that the cusp is preserved over a few or- 
bital periods, there are also claims that it can be flattened 
(Hayashi et al. 2003). One may also question the assumption 
of a cuspy initial profile in the light of recent refinements in 
studies of density profiles of dark matter haloes which find 
a flattening profile in the centre (Navarro et al. 2004; Stoehr 
2005). We believe however that, from the point of view of 
the present study, the inner slope of the halo is not really 
important since it should not affect our conclusions concern- 
ing the total mass or anisotropy. In the study of the dark 
matter distribution in the Coma cluster Lokas & Mamon 
(2003) considered a generalized profile with an arbitrary in- 
ner slope and found almost the same best-fitting mass and 
anisotropy for flat and cuspy inner profiles. 



5 RESULTS AND DISCUSSION 

In this Section, we report our attempts to reproduce the 
observed velocity moments shown in Fig. [5] using models 
described by equations with the mass distribution 

given by the sum of the two contributions 1201 and 1221 
discussed in the previous Section: 



M(r) = M,(r) + M A {r) 



(23) 



The density profile v{r) of the tracer population of stars 
is given by equation 11911 and the surface brightness I(R) 
by the Sersic formula HI 711 with m = 0.83. As discussed 
by Lokas & Mamon (2003) and Sanchis et al. (2004), while 
studying velocity dispersion is useful to constrain the mass, 
the kurtosis is mostly sensitive to the velocity anisotropy. 
However, as will become clear below, when the anisotropy is 
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constrained, the degeneracy between the mass distribution 
and velocity anisotropy, usually present when only velocity 
dispersion is considered, is broken and the parameters de- 
scribing the mass distribution can also be constrained. 

We begin by fitting the line-of-sight velocity dispersion 
profiles shown by the upper curves in Fig. [2] and adjusting 
three parameters: the anisotropy parameter j3, total dark 
mass in units of the total stellar mass Md/Ms and the break 
radius of the dark matter profile in units of the Sersic radius 
rb/Rs measuring the extent of the dark matter halo. The 
results for our different samples of 189, 203 and 207 stars 
are shown in Fig. [3] so that the best-fitting Mrj /Ms (up- 
per panel) and the velocity anisotropy parameter /3 (middle 
panel) are displayed as a function of r^/Rs- The thinner 
lines in the lower panel of Fig. [3] show the corresponding 
goodness of fit \ 2 ■ The lower panel of Fig. [3] proves that 
none of the parameters can be constrained from the analy- 
sis of velocity dispersion alone: \ 2 ls na t for a large range of 
»~b /Rs and does not discriminate between different values of 
the parameters. 

As discussed in Section 3, the sampling distributions of 
<ti os and (logKios) 1 / 10 are independent to a good approxima- 
tion, hence we can use the same \ 2 minimization scheme to 
find joint constraints following from fitting both quantities. 
Before that, however, we calculate \ 2 using the total of 10 
data points for each sample of 189, 203 and 207 stars for 
the values of the parameters dictated by the best fit to the 
velocity dispersion data in Fig. |3] The results are shown in 
the lower panel of the Figure with thicker lines. We imme- 
diately see that \ 2 has a minimum only for our sample with 
189 stars. For the other samples its values decrease for in- 
creasing rb / Rs and Mb /Ms and the minimum is not reached 
even if we go to incredibly high masses like Mo /Ms ~ 10 6 . 
Such masses would be comparable to the mass of the Milky 
Way itself. 

The reason for this behaviour can be understood by 
referring to Fig. 21 where we plot the predicted velocity 
moments for different values of the anisotropy parameter 
p. For simplicity we assumed the mass distribution as ob- 
tained from the fit to the sample of 189 stars (see eq. 1241 
below). It should be kept in mind that the kurtosis profile 
depends mainly on anisotropy, while the velocity dispersion 
is degenerate in anisotropy versus mass (taking a more ex- 
tended mass distribution for the plot would result in a more 
strongly increasing <ti os in analogy to more tangential orbits 
but the kurtosis profile would remain unaltered). As seen in 
the Figure and already discussed by Lokas & Mamon (2003) 
and Sanchis et al. (2004) (see also Gerhard 1993), tangen- 
tial orbits produce a decreasing kurtosis profile, as is indeed 
observed for the sample with 189 stars. But the kurtosis val- 
ues of samples with 207 and 203 stars are rather high (close 
to Gaussian) and therefore f3 values close to isotropic or 
mildly radial are preferred while increasing velocity disper- 
sion profiles require tangential orbits (see the middle panel 
of Fig. Without tangential orbits such velocity disper- 
sion profiles could only be reproduced if the mass profiles 
are extended and this is the reason the x 2 values decrease 
for higher masses and break radii. This indicates that the 
samples with 207 and 203 stars should be treated with cau- 
tion. In the case of the sample with 203 stars a reasonably 
good fit to both moments can be obtained, although for 
rather high masses. For the sample with 207 stars however, 
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Figure 3. Results of fitting only the line-of-sight velocity disper- 
sion data. The best fitting parameters Md/Ms and /3 are shown 
in the two upper panels as a function of rj, /Rs for different stellar 
samples. The thinner lines in the lower panel give the goodness 
of fit x 2 from fitting velocity dispersion. The thicker lines show 
X 2 values obtained for the same parameters if the kurtosis data 
are included. 
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Figure 4. Linc-of-sight velocity dispersion <ri oa (upper panel) and 
kurtosis (lower panel) as a function of projected radius, for differ- 
ent anisotropies: radial (/3 = 1), isotropic (/3 = 0) and moderately 
tangential (/3 = —4, i.e. crg/a r = 2.2). The mass distribution is 
given by the best-fit for the sample with 189 stars, eq. 1241 . 
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Figure 5. Line-of-sight velocity dispersion cri os (upper panel) 
and dimensionless line-of-sight kurtosis parameter (log Kios) 1 ^ 10 
(lower panel) with la error bars for the sample with 189 stars. 
The curves represent the best-fitting model. The thin horizon- 
tal line in the lower panel marks the value of (log3) 1//10 = 0.93 
corresponding to the Gaussian velocity distribution. 



the fit is bad and improves only for incredibly high masses, 
which may indicate that this sample yields inconsistent ve- 
locity dispersion and kurtosis profiles and therefore that it 
almost certainly includes unbound stars. 

One may ask to what extent this conclusion depends on 
our assumed constant anisotropy. It has recently been veri- 
fied using the cosmological iV-body simulations (Wojtak et 
al. 2005) that the radial velocity moments of bound objects, 
eqs. (|1J and ||KJ, are actually sensitive to the local value of 
the anisotropy parameter. We expect this to be true also for 
the projected moments since the dominant contribution to 
the projected moment near the centre of the object comes 
from the stars that are actually near the centre because they 
are more numerous than the foreground and background 
stars, while the value of the projected moment at some pro- 
jected distance R from the centre is affected only by stars 
located at true distances r > R. This behaviour in the case 
of velocity dispersion has been explicitly verified by Lokas & 



Mamon (2001) for the Osipkov-Merritt anisotropy (see their 
Figures 1 and 6). We therefore expect that the anisotropy 
parameter changing e.g. from isotropic orbits near the cen- 
tre to tangential orbits at larger distances would produce 
cTi os and Kios profiles following the dashed line at small R 
and solid line at large R in the upper and lower panel of 
Fig. ^respectively. The kurtosis profile would then decrease 
even more strongly and the conclusion concerning the sam- 
ples with 203 and 207 stars would be similar as in the case 
of constant (5. Interestingly, the dispersion profile interpo- 
lating between the dashed and solid line in the upper panel 
of Fig. [I] could reproduce the dip in the velocity dispersion 
profile seen in the data for all samples (see Fig. |5J and pre- 
served for different binnings of the stars. 

We therefore conclude that the choice of the sample of 
stars is the dominant source of uncertainty when trying to 
determine the mass distribution in Draco. In the following 
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we will only consider the sample with 189 stars and investi- 
gate the next important source of uncertainty following from 
the sampling errors we assigned to our data points. The er- 
rors in the adopted parameters of Draco concerning the sur- 
face distribution of stars, luminosity or distance probably 
have much smaller effect (see Lokas 2002). 

The joint fitting of the velocity dispersion and kurtosis 
for the sample with 189 stars gives results not very different 
from those indicated by the minimum of y 2 (thicker dotted 
line) in the lowest panel of Fig. |3] The best-fitting parame- 
ters (with x 2 /N = 7.4/10) are 

Mn/Ms = 99, r h /R s = 1.9, (3 = -0.4. (24) 

The best estimate of the total mass of the Draco dwarf is 
then 6.6 x 10 7 Mq . We emphasize that the result is not robust 
and could be very different if different stellar samples were 
taken. Our choice of velocities has been rather arbitrary and 
one could probably get reasonable results also with slightly 
larger or smaller samples. One could for example construct 
samples by removing the stars with most discrepant veloc- 
ities one by one and performing a similar fitting procedure 
as described here. One would then get a series of best fitting 
masses of decreasing values. 

The mass estimate found above is significantly smaller 
than the ones obtained with the NFW mass distribution 
and older data (Lokas 2002), which were even of the order 
of 10 9 Mq for reasonable concentrations. This was due to 
the more extended nature of the NFW profile, which was 
modified here by adding an exponential cut-off. Recently, it 
has been proposed that, if the dwarf spheroidals of the Local 
Group indeed possess such high masses, then their low abun- 
dance could be understood by referring to the cosmological 
mass function, which predicts many fewer objects with such 
high mass (Hayashi et al. 2003). Our present lower mass esti- 
mate could be consistent with this solution to the overabun- 
dance problem if the dwarfs originally indeed had such high 
masses but were then tidally stripped in the Milky Way po- 
tential. Such a scenario could also explain why dwarfs man- 
aged to build up their stellar content in presently shallow 
potential wells (Kravtsov, Gnedin & Klypin 2004). 

Our estimate of Draco's mass is not directly compara- 
ble to that of Kleyna et al. (2002) since they used differ- 
ent dark matter models and a smaller stellar sample (con- 
structed from an older data set but in a weakly restrictive 
way similar to our sample with 207 stars). Their best es- 
timate of the mass inside three core radii was 8 x 10 7 Mq, 
which is, as would be expected, already more than our total 
mass value. 

The velocity moments obtained with the sets of param- 
eters listed in (1241 are shown in Fig. together with the 
data for 189 stars. The dark matter density profile following 
from equation 12111 with those parameters is plotted in Fig.||j] 
together with the stellar mass density profile p*(r) = Tu(r) 
with v(r) given by equation 1191 . Fig. |S| also shows the cu- 
mulative mass-to-light ratio, i.e. the ratio of the total mass 
distribution to the luminosity distribution in the V-band 
inside radius r 

M/Ly = ¥£L, (25) 
L v (r) 

where M(r) is the sum of two contributions from stars and 
dark matter given by equation 1231 with 12QH and I22H , while 




r [arcmin] 

Figure 6. Best fitting dark matter density profile (solid line) in 
comparison to the stellar mass density profile (dashed line). The 
dotted line shows the cumulative mass-to-light ratio for the best 
fitting profile. 

Ly(r) = M*(r)/Tv- The behaviour of this quantity is easily 
understood: the dark matter cusp is steeper than that of the 
luminosity distribution therefore the ratio increases towards 
the centre. This is due to our adopted profile of dark matter 
halo, which, according to the simulations of Kazantzidis et 
al. (2004), preserves the NFW cusp. With a flatter dark 
matter inner slope (as recently proposed by Navarro et al. 
2004 and Stoehr 2005 for global halos, and by Hayashi et 
al. 2003 for tidally stripped ones), the data may well be 
reconciled with a constant mass-to-light ratio. 

At large distances, the mass-to-light ratio converges to 
a constant value Tv(l + Md/Ms) = 299 since both lu- 
minosity and dark matter distributions fade exponentially. 
This result is comparable to the estimate by Kleyna et al. 
(2002) who find M/Lv — 330 within three core radii, while 
Odenkirchen et al. (2001), assuming that mass follows light, 
obtain M/Li = 146. Our result is also consistent with the 
study of Lokas (2002) using older data, where dark matter 
was modelled with a generalized NFW distribution. There 
M/Ly > 100 was found at R > 10 arcmin, and the mass-to- 
light ratio was increasing at large radii due to the assump- 
tion of an untruncated, p oc r -3 , dark matter envelope. 

Fig. [7] illustrates the uncertainties in the estimated pa- 
rameters due to the sampling errors of the velocity moments. 
Here, we plot the cuts through the confidence region in three 
possible parameter planes with probability contours corre- 
sponding to la (68 percent confidence), 2a (95 percent) and 
3a (99.7 percent) i.e. Ax 2 = X 2 ~ Xmin = 3-53, 8.02, 14.2, 
where Xmin — 7.4. Contrary to the case where only the ve- 
locity dispersion was fitted, some constraints on parameters 
can now be obtained, but the allowed ranges of parameters 
are large. This is especially the case for parameters Mo /Ms 
and rb/Rs which are highly degenerate. Since the parame- 
ters Md/Ms and r^/Rs are very weakly correlated with /?, 
as can be seen from the two upper panels of Fig. |7| the cuts 
provide a good approximation of the true uncertainties in 
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the parameters. However, these errors, due to the low num- 
ber of velocities per bin, are much smaller than uncertainties 
due to our choice of stars to include in the sample. 

As already stated above, we find the preferred stellar 
orbits to be weakly tangential. For the sample with 189 
stars we get the best-fitting = —0.4 (corresponding to 
oe/oy « 1.2) and the expected uncertainties can be read 
from the two upper panels of Fig. Q Note that according 
to the Figure the orbits are consistent with isotropy at 1 <r 
level, while a strongly radial or a strongly tangential ve- 
locity distribution are ruled out. It must be remembered, 
however, that these conclusions were reached with the as- 
sumption of negligible binary fraction in Draco. As discussed 
in Section 3, the correction for binaries significantly affects 
the values of kurtosis (see Fig. but again more so for 
the sample with 189 stars than for the other samples which 
have higher values of kurtosis. If the binary fraction in Draco 
is significant, e.g. of the order of 7 = 0.5 as for the solar 
neighbourhood according to current estimates (Duquennoy 
& Mayor 1991; De Rijcke & Dejonghe 2002) then the intrin- 
sic kurtosis will be decreased and the inferred orbits even 
more tangential. Reproducing the velocity dispersion profile 
(which is not strongly affected by the binaries) would then 
require a less extended mass distribution and the inferred 
masses will be even smaller. 

We have mentioned before that the main gain from in- 
cluding the kurtosis in the analysis is the ability to con- 
strain the anisotropy of stellar orbits. Interestingly, weakly 
tangential orbits in Draco appear to be a robust conclusion 
from our study; whatever the sample of the stars taken and 
binary fraction assumed the preferred models have f3 < 
or og/a r > 1. This result seems to contradict the general 
belief, based e.g. on iV-body simulations, that virialized sys- 
tems should have orbits close to isotropic or even mildly 
radial (see Mamon & Lokas 2005). However, as discussed 
by Kazantzidis et al. (2004), tidally stirred stellar popula- 
tions can develop tangential orbits even when evolving in an 
isotropic halo. 

A possible scenario for the formation of dwarf 
spheroidal galaxies has been proposed by Mayer et al. 
(2001) who used high-resolution iV-body/SPH simulations 
to study the evolution of rotationally supported dwarf ir- 
regular galaxies moving on bound orbits in the massive 
dark matter halo of the Milky Way. They showed that the 
tidal field induces severe mass loss in their halos and the 
low surface brightness dwarfs are transformed into pressure- 
supported dwarf spheroidals with preferably tangential or- 
bits outside the centre (see Fig. 25 of the preprint version of 
Mayer et al. 2001). 

Besides, tangential orbits may not be restricted to dwarf 
spheroidals. A recent study of kinematics of a sample of 
Virgo Cluster dwarf ellipticals (Geha, Guhathakurta & van 
der Marel 2002) also shows the preference of tangential or- 
bits and actually the mechanism to produce them may be 
similar as in the case of dwarf spheroidals (Mayer et al. 
2001) only with high surface brightness galaxies as progen- 
itors. Also, elliptical galaxies formed by major mergers of 
spiral galaxies often show tangential anisotropy in their in- 
ner regions (Dekel et al. 2005), which is caused by some of 
the stars having formed within the gaseous disk, which is in 
nearly circular rotation in the merger remnant. 
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Figure 7. Cuts through the 1 a, 2cr and 3 <r probability con- 
tours in the parameter space obtained from fitting <ti os and 
(log Kiog) 1 / 10 for the sample with 189 stars. Dots indicate the 
best-fitting parameters. 
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Figure 8. Expected contamination of Draco by Milky Way stars. 
The diagram shows the expected total number of Milky Way stars 
of given velocities (higher histogram, thin solid line) and the cor- 
responding number if we restrict ourselves to the red giant branch 
of the color-magnitude diagram of Draco (lower histogram, thick 
solid line). The dashed line shows the velocity distribution of 
stars actually observed in the direction of Draco by Wilkinson 
et al. (2004). The contamination by giants in the velocity range 
of Draco, indicated by two vertical dashed lines, is very small. 



6 INTERLOPERS FROM THE MILKY WAY 
OR TIDAL DEBRIS? 

We have shown that the main issue in the analysis of the 
velocity moments is the question of which stars should be in- 
cluded in the analysis, or, to put it in a different way, which 
of them are unbound and only artificially affect the mo- 
ments, e.g. inflate the velocity dispersion. We argued that 
the presence of unbound stars leads to inconsistent veloc- 
ity dispersion and kurtosis which cannot both be fitted by 
equilibrium models following from the Jeans formalism. It 
is therefore essential to study the possible origin of those 
unbound stars. 

An immediate first guess would be that the stars are 
the members of the Milky Way. In order to verify this hy- 
pothesis we have ran the Besancon Galaxy model (Robin et 
al. 2003, http://bison.obs-besancon.fr/modele/ I in the di- 
rection of Draco, which gives the velocity distribution shown 
by the higher histogram (thinner solid line) in Fig. |H| If stars 
of all types are considered, as they are for this histogram, 
then the expected number of Milky Way interlopers with 
Draco-like velocities in the range -325 km s _1 < v < -262 
km s _1 , indicated by the two vertical dashed lines in the Fig- 
ure, would be 15 ± 0.25. However, if we restrict the analysis 
to the red giant branch (RGB) of the Draco color-magnitude 
diagram with ±0.075 mag from the typical trend we end up 
with a much lower expected number density of stars repre- 
sented by the lower histogram (thicker solid line) in Fig. |H| 
Since the observations of Draco stars by Kleyna et al. (2002) 
and Wilkinson et al. (2004) were indeed restricted to RGB 
stars we expect the contamination from the Milky Way stars 
to be very small with 1.42 ± 0.17 expected interlopers. The 



velocity distribution of all RGB stars observed by Wilkin- 
son et al. (2004) is shown with a dashed line in the plot. 
We can see that the level of this histogram right outside 
the range of velocities adopted for Draco stars is somewhat 
higher than expected from the Besancon model. A naive 
interpolation of the observed velocity distribution (dashed 
histogram in Fig- |HJ outside of the Draco velocity zone would 
yield 6.5 ± 2.3 interlopers from the Milky Way. 

Given our finding that the kinematical modelling of the 
189-star sample yields consistent results, while those of the 
207- and 203-star samples do not, the total number of inter- 
lopers in Draco is of order of 207— 189 = 18 stars. Therefore, 
given our estimates of the previous paragraph of the num- 
ber of Milky Way interlopers, an important fraction (prob- 
ably the majority) of unbound stars in the neighbourhood 
of Draco dwarf must originate in Draco itself. An obvious 
mechanism for the production of such stars is the tidal strip- 
ping of Draco by the Milky Way potential. Although no clear 
signature of tidal tails has been identified in the luminosity 
distribution of Draco (e.g. Klessen et al. 2003), it is still 
quite possible that the tails are oriented along the line of 
sight. 

Wilkinson et al. (2004) have recently argued that the 
tidal scenario is only feasible if the mass of Draco is below 
5 x 10 7 Mq and its dark matter distribution is low outside 
roughly 30 arcmin. Interestingly, the results of our analysis 
seem to almost fulfill these requirements: our best estimate 
of the mass for the sample with 189 stars is 6.6 x 10 7 Mq and 
the dark matter distribution indeed fades at larger distances, 
as demonstrated by Fig.|S| Wilkinson et al. (2004) also point 
out that such a tidally detached population of stars should 
be heated up, which is contradicted by the decreasing veloc- 
ity dispersion of Draco at large distances, unless only cold 
stars remain which are then recaptured by the dwarf. We 
find that the fall of the velocity dispersion at large angular 
distances from the centre of Draco is not a very strong fea- 
ture and its visibility depends on the chosen radial binning 
of stars. 



7 CONCLUDING REMARKS 

Given the new data on the Draco dwarf spheroidal, espe- 
cially the velocity measurements by Kleyna et al. (2002) and 
Wilkinson et al. (2004), we have attempted to constrain the 
dark matter distribution in Draco by modelling not only the 
dispersion, but also the kurtosis of the line-of-sight veloc- 
ity distribution. The analysis of both moments allows us to 
break the degeneracy between the mass distribution and ve- 
locity anisotropy usually present in the analyses of velocity 
dispersion. Still, many uncertainties remain in the analysis. 
We show that the results can be very different depending 
on the sample of stars chosen, i.e. which stars are included 
in the sample and which are treated as interlopers. Besides, 
due to the limited number of measured velocities, the sam- 
pling errors of the moments are large and the parameters 
cannot be strongly constrained. 

There are few possible remedies to improve the present 
situation: to have more velocity measurements, to further 
constrain the dark matter profiles on theoretical grounds 
and to model the influence of interlopers in more detail. 
All these approaches seem feasible in the near future. On 
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the observational side, the measurements of a few hun- 
dred velocities per dwarf galaxy are planned or already be- 
ing performed using large telescopes with multi-fiber spec- 
troscopy like Magellan (M. Mateo, private communication). 
On the theoretical side, progress might be even faster with 
the rapidly increasing resolution of the iV-body simulations. 
The recent result of Kazantzidis et al. (2004) itself which 
we have used in the present work relied on increased reso- 
lution to show that the cusp of the dark matter profile is 
not flattened due to tidal interaction as was claimed earlier 
(Stoehr et al. 2002; Hayashi et al. 2003). The result however 
depended to some extent on the parameters of the orbit of 
the dwarf around the Milky Way, which in the case of Draco 
are poorly known. 

As for the third issue, the treatment of interlopers or 
unbound stars, in our opinion this is at present the most 
important problem in the analysis of the dark matter dis- 
tribution in dwarf spheroidal galaxies. We have shown that 
the supposed interlopers can alter the profiles of velocity 
moments and significantly change the estimated parameters 
of the dark matter distribution. 
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